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Abstract 

A higher order Godunov method for the radiation subsystem of radiation hydrody- 
namics is presented. A key ingredient of the method is the direct coupling of stiff 
source term effects to the hyperbolic structure of the system of conservation laws; it 
is composed of a predictor step that is based on Duhamel's principle and a corrector 
step that is based on Picard iteration. The method is second order accurate in both 
time and space, unsplit, asymptotically preserving, and uniformly well behaved from 
the photon free streaming (hyperbolic) limit through the weak equilibrium diffusion 
(parabolic) limit and to the strong equilibrium diffusion (hyperbolic) limit. Numerical 
tests demonstrate second order convergence across various parameter regimes. 



1 Introduction 

Radiation hydrodynamics is a fluid description of matter (plasma) that absorbs and emits 
electromagnetic radiation and in so doing modifies dynamical behavior. The coupling be- 
tween matter and radiation is significant in many phenomena related to astrophysics and 
plasma physics, where radiation comprises a major fraction of the internal energy and mo- 
mentum and provides the dominant transport mechanism. Radiation hydrodynamics governs 
the physics of radiation driven outflows, supernovae, accretion disks, and inertial confinement 
fusion [UI2]- Such physics is described mathematically by a nonlinear system of conservation 
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laws that is obtained by taking moments of the Boltzmann and photon transport equations. 
A key difficulty is choosing the frame of reference in which to take the moments of the 
photon transport equation. In the comoving and mixed frame approaches, one captures the 
matter /radiation coupling by adding relativistic source terms correct to 0(u/c) to the right- 
hand side of the conservation laws, where u is the material flow speed and c is the speed of 
light. These source terms are stiff because of the variation in time/length scales associated 
with such problems [3j. This stiffness causes numerical difficulties and makes conventional 
methods such as operator splitting and method of lines breakdown jH [5] . 

Previous research in numerically solving radiation hydrodynamical problems was carried 
out by Castor 1972, Pomraning 1973, Mihalas & Klein 1982, and Mihalas & Mihalas 1984 
[2, El E], [7] . There are a variety of algorithms for radiation hydrodynamics. One of the sim- 
plest approaches was developed by Stone, Mihalas, & Norman 1992 and implemented in the 
ZEUS code, which was based on operator splitting and Crank-Nicholson finite differencing 
[8]. Since then, higher order Godunov methods have emerged as a valuable technique for 
solving hyperbolic conservation laws (e.g., hydrodynamics), particularly when shock captur- 
ing and adaptive mesh refinement is important [9]. However, developing upwind differencing 
methods for radiation hydrodynamics is a difficult mathematical and computational task. 
In many cases, Godunov methods for radiation hydrodynamics either: (i) neglect the het- 
erogeneity of weak/strong coupling and solve the system of equations in an extreme limit 
[TOl [TTj . (ii) are based on a manufactured limit and solve a new system of equations that 
attempts to model the full system [T2l [13] , or (Hi) uses a variation on flux limited diffusion 
[T4"t [15]. All of these approaches do not treat the full generality of the problem. For ex- 
ample, in a series of papers, Balsara 1999 proposed a Riemann solver for the full system of 
equations [TE] . However, as pointed out by Lowrie & Morel 2001, Balsara's method failed to 
maintain coupling between radiation and matter. Moreover, Lowrie & Morel were critical of 
the likelihood of developing a Godunov method for full radiation hydrodynamics [IT] . 

In radiation hydrodynamics, there are three important dynamical scales and each scale is 
associated with either the material flow (speed of sound), radiation flow (speed of light), 
or source terms. When the matter-radiation coupling is strong, the source terms define the 
fastest scale. However, when the matter-radiation coupling is weak, the source terms define 
the slowest scale. Given such variation, one aims for a scheme that treats the stiff source 
terms implicitly. Following work by Miniati & Colella 2007, this paper presents a method 
that is a higher order modified Godunov scheme that directly couples stiff source term effects 
to the hyperbolic structure of the system of conservation laws; it is composed of a predic- 
tor step that is based on Duhamel's principle and a corrector step that is based on Picard 
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iteration [IB] . The method is explicit on the fastest hyperbolic scale (radiation flow) but 
is unsplit and fully couples matter and radiation with no approximation made to the full 
system of equations for radiation hydrodynamics. 

A challenge for the modified Godunov method is its use of explicit time differencing when 
there is a large range in the time scales associated with the problem, c/a^ ^> 1 where 
is the reference material sound speed. One could have built a fully implicit method 
that advanced time according to the material flow scale, but a fully implicit approach was 
not pursued because such methods often have difficulties associated with conditioning, are 
expensive because of matrix manipulation and inversion, and are usually built into central 
difference schemes rather than higher order Godunov methods. An explicit method may even 
out perform an implicit method if one considers applications that have flows where cja^ < 
10. A modified Godunov method that is explicit on the fastest hyperbolic scale (radiation 
flow) as well as a hybrid method that incorporates a backward Euler upwinding scheme for 
the radiation components and the modified Godunov scheme for the material components 
are under construction for full radiation hydrodynamics. A goal of future research is to 
directly compare these two methods in various limits for different values of c/aoo. 



2 Radiation Hydrodynamics 

The full system of equations for radiation hydrodynamics in the Eulerian frame that is 
correct to 0(1/ C) is: 
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P r = \E r (closure relation). 
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For the material quantities, p is density, m is momentum, p is pressure, E is total energy den- 
sity, and T is temperature. For the radiative quantities, E r is energy density, F r is flux, P r is 
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pressure, and f is the variable tensor Eddington factor. In the source terms, o a is the absorp- 
tion cross section, a s is the scattering cross section, and a t = a a + a s is the total cross section. 

Following the presentation of Lowrie, Morel, & Hittinger 1999 and Lowrie & Morel 2001, the 
above system of equations has been non-dimensionalized with respect to the material flow 
scale so that one can compare hydrodynamical and radiative effects as well as identify terms 
that are O(u/o). This scaling gives two important parameters: C = c/aoo, P = . C 
measures relativistic effects while P measures how radiation affects material dynamics and 
is proportional to the equilibrium radiation pressure over material pressure. a r = r^Jh 
is a radiation constant, is the reference material temperature, and is the reference 
material density. 

For this system of equations, one has assumed that scattering is isotropic and coherent in 
the comoving frame, emission is defined by local thermodynamic equilibrium (LTE), and 
that spectral averages for the cross-sections can be employed (gray approximation). The 
coupling source terms are given by the modified Mihalas-Klein description [TT], [19] which is 
more general and more accurate than the original Mihalas-Klein source terms [3] because it 
maintains an important 0(1/C 2 ) term that ensures the correct equilibrium state and relax- 
ation rate to equilibrium [TTJ [19] . 

Before investigating full radiation hydrodynamics, it is useful to examine the radiation sub- 
system, which is a simpler system that minimizes complexity while maintaining the rich 
hyperbolic-parabolic behavior associated with the stiff source term conservation laws. This 
simpler system allows one to develop a reliable and robust numerical method. Consider 
Equations HI [5] for radiation hydrodynamics in one spatial dimension not affected by trans- 
verse flow. If one only considers radiative effects and holds the material flow stationary 
such that it — > 0, then the conservative variables, fluxes, and source terms for the radiation 
subsystem are given by: 

^ + C§=C,.<I--ft), (7) 

Motivated by the asymptotic analysis of Lowrie, Morel, & Hittinger 1999 for full radiation 
hydrodynamics, one investigates the limiting behavior for this simpler system of equations. 
For non- relativistic flows 1/C = 0(e), where e C 1. Assume that there is a moderate amount 
of radiation in the flow such that P = 0(1). Furthermore, assume that scattering effects are 
small such that o s jo t = 0(e). Lastly, assume that the optical depth can be represented as 
£ = i'mat/At = £ mat <7 t , where \ t is the total mean free path of the photos and £ mat = 0(1) 
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is the material flow length scale [T9j . 

Free Streaming Limit a a ,a t ~ 0(e): In this regime, the right-hand-side of Equations [7] 
and M is negligible such that the system is strictly hyperbolic. / — > 1 and the Jacobian of 
the quasilinear conservation law has eigenvalues ±C: 



Weak Equilibrium Diffusion Limit cx a ,cr t ~ 0(1) : One obtains this limit by plugging 
in a a ,a t ~ 0(1), matching terms of like order, and combining the resulting equations. From 
the definition of the equilibrium state, E r = T 4 and F r = Therefore, the system is 

parabolic and resembles a diffusion equation, where / — > 1/3: 



dE r _ C d 2 E r 
dt 3cr t dx 2 



1 8E r 

= ~3F t -W- < 12 > 

Strong Equilibrium Diffusion Limit a a , o t ~ 0(1/ e) : One obtains this limit by plugging 
in cr a ,<7f ~ 0(l/e) and following the steps outlined for the weak equilibrium diffusion limit. 
One can consider the system to be hyperbolic, where / — > 1/3 and the Jacobian of the 
quasilinear conservation law has eigenvalues ±e: 

dE r , . 

ir= - (13 > 

F r = 0. (14) 

Lowrie, Morel, & Hittinger 1999 investigated an additional limit for full radiation hydro- 
dynamics, the isothermal regime. This limit has some dynamical properties in common 
with the weak equilibrium diffusion limit, but its defining characteristic is that the material 
temperature T(x,t) is constant. When considering the radiation subsystem, there is little 
difference between the weak equilibrium diffusion and isothermal limits because the mate- 
rial quantities, including the material temperature T, do not evolve. T enters the radiation 
subsystem as a parameter rather than a dynamical quantity. 

3 Higher Order Godunov Method 

In one spatial dimension, systems of conservation laws with source terms have the form: 

dU , dF(U) 

-bT + ^x- = Siu) > (15) 
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where U : R x R — > R n is an n-dimensional vector of conserved quantities. For the radiation 
subsystem: 

\F r )' \ CfE r )' \CS F ) y -Co,F r 

The quasilinear form of this system of conservation laws is: 

dU A dU 0/rn , dF ( C \ 

m +A lT x = s ^ A =du=[cf o J" (16) 

A has eigenvalues A = ±/ 1//2 C as well as right eigenvectors R (stored as columns) and left 
eigenvectors L (stored as rows): 



R 




/ \ 1/2 \ 

1 _ 1 / 1 \ ' \ 

2 2 \f 

, ' v'l/ 2 

1 1^1 

2 2 



(17) 



V J 



Godunov's method obtains solutions to systems of conservation laws by using characteristic 
information within the framework of a conservative method: 

U? +1 = U?~^ (Fi + x,2 ~ F t .- 1/2 ) + AtS(Un- (18) 

Numerical fluxes ^±1/2 are obtained by solving the Riemann problem at the cell interfaces 
with left/right states to get U^j^ and computing Fi±i/ 2 = ^(11^^)^ where i represents 
the location of a cell center, i ±1/2 represents the location cell faces to the right and left of 
i, and superscripts represent the time discretization. An HLLE (used in this work) or any 
other approximate Riemann solver may be employed because the Jacobian dF/dU for the 
radiation subsystem is a constant valued matrix and by definition a Roe matrix pfl, El 120] . 
This property also implies that one does not need to transform the system into primitive 
variables (VyW). The power of the method presented in this paper is that the spatial 
reconstruction, eigen- analysis, and cell-centered updating directly plug into conventional 
Godunov machinery. 

3.1 Predictor Step 

One computes the flux divergence (V • _F) n+1 / 2 by using the quasilinear form of the system 
of conservation laws and the evolution along Lagrangian trajectories: 

DU + = s{v)i a l = A _ uh EE = au + /a\ K 



Dt dx v n ' Dt dt V dx. 

From the quasilinear form, one derives a system that includes (at least locally in time and 
state space) the effects of the stiff source terms on the hyperbolic structure. Following the 
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analysis of Miniati & Colella 2007 and Trebatich et al 2005 [T8l |2T] . one applies Duhamel's 
principle to the system of conservation laws, thus giving: 



DU 



cff 



Dt 



T sM [~A L ^ + S?), (20) 



where X s is a propagation operator that projects the dynamics of the stiff source terms 
onto the hyperbolic structure and S n = VuS\u n - The subscript n designates time t = t n . 
Since one is considering a first order accurate predictor step in a second order accurate 
predictor-corrector method, one chooses rj = At/2 and the effective conservation law is: 

^+X Sn (At/2)A L ^ = l Sn (At/2)S n , =► ^ + A cS ^- = l Sn (At/2)S n , (21) 

where A e s = X Sn (At/2)A L + ul. In order to compute X s , one first computes S n . Since C, 
cx a , and a t are constant and one assumes that J^, J^- = 0: 







Ig is derived from Duhamel's principle and is given by: 

1 /"At/2 

T Sn (At/2) = — I e^dr (23) 

1 _ e -«> a At/2 J _ e -C<7 t At/2 

a = £a a At/2 ' ^ = ' (24) 

Before applying Xg to A^, it is important to understand that moving-mesh methods can 
be accommodated in non-relativistic descriptions of radiation hydrodynamics whenever an 
Eulerian frame treatment is employed. These methods do not require transformation to the 
comoving frame [T7]. Since the non-dimensionalization is associated with the hydrodynamic 
scale, one can use M mes h = u from Lagrangean hydrodynamic methods. 

The effects of the stiff source terms on the hyperbolic structure are accounted for by trans- 
forming to a moving-mesh (Lagrangean) frame Al = A — ul, applying the propagation 
operator Xg to Al, and transforming back to an Eulerian frame A e g = Xg Al + ul [18J. 
However, because only the radiation subsystem of radiation hydrodynamics is considered 
«mesh = u — > 0. Therefore, the effective Jacobian is given by: 

^=i,; c ; c i, (25) 
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which has eigenvalues A e fr = ±(a/?) 1 / 2 / 1 / 2 C with the following limits: 

a a-, a t —> a,f3 —y 1 =>- A c ff — > ±/ 1 / 2 C, (free streaming) 
CajCt ^ oo =>- a, /3— >0 =>- A c flf — > ±e, (strong equilibrium diffusion). 

A e ff has right eigenvectors i? e ff (stored as columns) and left eigenvectors L e s (stored as rows) : 



/ i _if«V /a \ 

2 2 U/ 



(?) ' (f) 



I 1 ( °L 
\ 2 2^/3/ 



1/2 



(26) 



3.2 Corrector Step 



The time discretization for the source term is a single-step, second order accurate scheme 
based on the ideas from Dutt et al 2000, Minion 2003, and Miniati & Colella 2007 [TBI [22| [23] . 
Given the system of conservation laws, one aims for a scheme that has an explicit approach 
for the conservative flux divergence term V ■ F and an implicit approach for the stiff source 
term S(U). Therefore, one solves a following collection of ordinary differential equations at 
each grid point: 

= S(U) - (V ■ F) n+1 '\ (27) 

at 

where the time-centered flux divergence term is taken to be a constant source which is 
obtained from the predictor step. Assuming time t = t n , the initial guess for the solution at 
the next time step is: 

U = U n + At(I - AtV C /S(f/)| c/ n)- 1 (S(f/ n ) - (V • F) n+1 / 2 ), (28) 

where: 

(/ - AtVuSiU))- 1 = ( 1+A ' C ^ ° j. (30) 

The error e is defined as the difference between the initial guess and the solution obtained 
from the Picard iteration equation where the initial guess was used as a starting value: 

e(At) = U n + -y (S(U) + S(U n )^j - At(V • F) n+1 ' 2 - U. (31) 
Following Miniati & Colella 2007, the correction to the initial guess is given by [18J: 

6 (At) = (/ - AtVuSiU)^)- 1 e(At). (32) 
Therefore, the solution at time t = t n + At is: 

U n+l = U + 5(At). (33) 
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3.3 Stability and Algorithmic Issues 

The higher order Godunov method satisfies important conditions that are required for 
numerical stability [TS]. First, A c g = ±(a/?) 1//2 / 1 / 2 C indicates that the sub characteris- 
tic condition for the characteristic speeds at equilibrium is always satisfied, such that: 
^ _ < Ks < A < X^ s < A + . This condition is necessary for the stability of the system 
and guarantees that the numerical solution tends to the solution of the equilibrium equation 
as the relaxation time tends to zero. Second, since the structure of the equations remains 
consistent with respect to classic Godunov methods, one expects the CFL condition to apply: 
max(|A*|)£i < 1, * = -, 0,+. 

Depending upon how one carries out the spatial reconstruction to solve the Riemann problem 
in Godunov's method, the solution is either first order accurate in space (piecewise constant 
reconstruction) or second order accurate in space (piecewise linear reconstruction). Piecewise 
linear reconstruction was employed in this paper, where left /right states (with respect to the 
cell center) are modified to account for the stiff source term effects [HI [24] : 

= U? + ^ X Sn {^ ^ S{Uf) + \{±I - P^AU,) (34) 

P ± (AU t ) = ( L eff • W) ■ R k cS - (35) 

±A fc >0 

Left /right one-sided slopes as well as cell center slopes are defined for each cell centered 
quantity t/j. A van Leer limiter is applied to these slopes to ensure monotonicity, thus giving 
the local slope AC/j. 



4 Numerical Tests 

Four numerical tests spanning a range of mathematical and physical behavior were carried 
out to gauge the temporal and spatial accuracy of the higher order Godunov method. The 
numerical solution is compared with the analytic solution where possible. Otherwise, a self- 
similar comparison is made. Using piecewise constant reconstruction for the left/right states, 
one can show that the Godunov method reduces to a consistent discretization in each of the 
limiting cases. 



The optical depth r is a useful quantity for classifying the limiting behavior of a system that 
is driven by radiation hydrodynamics: 

1CLX 

u t dx = <x t (x max - x min ), (36) 
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Optically thin/thick regimes are characterized by: 

t < 0(1) (optically thin) 
t > 0(1) (optically thick). 

In optically thin regimes (free streaming limit), radiation and hydrodynamics decouple such 
that the resulting dynamics resembles an advection process. In optically thick regimes 
(weak/strong equilibrium diffusion limit), radiation and hydrodynamics are strongly cou- 
pled and the resulting dynamics resembles a diffusion process. 



The following definitions for the n-norms and convergence rates are used throughout this 
paper. Given the numerical solution q r at resolution r and the analytic solution u, the error 
at a given point i is: e\ = q\ — u. Likewise, given the numerical solution q r at resolution r 
and the numerical solution q r+1 at the next finer resolution r + 1 (properly spatially averaged 
onto the coarser grid), the error resulting from this self-similar comparison at a given point 
i is: e\ = q\ — The 1-norm and max- norm of the error are: 

L 1 = y^\e r i \Ax r , L max = max|e-|. (37) 

i 

The convergence rate is measured using Richardson extrapolation: 

n ~ ln(Aa;7A^+ 1 ) ' 1 ' 

4.1 Exponential Growth/Decay to Thermal Equilibrium 

The first numerical test examines the temporal accuracy of how variables are updated in the 
corrector step. Given the radiation subsystem and the following initial conditions: 

= constant across space, F r ° = 0, T = constant across space, 

F r — > for all time. Therefore, the radiation subsystem reduces to the following ordinary 
differential equation: 

HE 

^ = Ca a (T*-E r ), (39) 
which has the following analytic solution: 

E r = T 4 + (E° r - T 4 )exp(-G7 a t). (40) 

For E® < T 4 and F r ° = 0, one expects exponential growth in E r until thermal equilibrium 
(En = T 4 ) is reached. For E® > T 4 and F r ° = 0, one expects exponential decay in E r 
until thermal equilibrium is reached. This numerical test allows one to examine the order of 
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accuracy of the stiff ODE integrator. 
Parameters: 

C = 10 5 , <t« = 1, a t = 2, / = 1, 
N^u = [32, 64, 128, 256], 
^min = 0, x max = 1, Ax = — , CFL = 0.5, At 

Wcell 

IC for Growth: £ r ° = 1, F r ° = 0, T = 10, 
1C for Decay: E° r = 10 4 , F r ° = 0, T = 1. 

1 

0.8 
* 0.6 
m ^0.4 
0.2 

°0 1 2 3 4 5 

t/o C 
a 

Figure 1: Exponential growth/decay to 
thermal equilibrium. N ce u = 256. 
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2.0 
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9.3E-3 


2.0 


9.3E-3 


2.0 


9.3E-3 


2.0 
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256 


2.3E-3 


2.0 


2.3E-3 


2.0 


2.3E-3 


2.0 


2.3E-3 


2.0 



Table 1: Errors and convergence rates for exponential growth/decay in E r to thermal 
equilibrium. Errors were obtained through analytic comparison, t = 10~ 5 = l/a a C 



From Figured], one sees that the numerical solution corresponds with the analytic solution. In 
Table 1, the errors and convergence rates are identical for growth and decay. This symmetry 
illustrates the robustness of the Godunov method. Furthermore, one finds that the method 
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is well behaved and obtains the correct solution with second order accuracy for stiff values 
of the e folding time ( 1 ^ t c > 1), although with a significantly larger amplitude in the norm 
of the error. This result credits the flexibility of the temporal integrator in the corrector step. 



In a similar test, the initial conditions for the radiation energy and flux are zero and the tem- 
perature is defined by some spatially varying profile (a Gaussian pulse). As time increases, 
the radiation energy grows into T(x) 4 . Unless the opacity is sufficiently high, the radiation 
energy approaches but does not equal T(x) 4 . This result shows that the solution has reached 
thermal equilibrium and any spatially varying temperature will diffuse. 



4.2 Free Streaming Limit 

In the free streaming limit, r <C 0(1) and the radiation subsystem reduces to Equations 
IU [KB If one takes an additional temporal and spatial partial derivative of the radiation 
subsystem in the free streaming limit and subtracts the resulting equations, then one finds 
two decoupled wave equations that have the following analytic solutions: 

E r (x,t) = E (x-f 1/2 Ct), (41) 
F r (x,t) =F (x-f 1/2 Ct). (42) 



Parameters: 



C = 10°, <r a = 10"°, a t = 10"°, / = 1, T = 1, 
N cM = [32, 64, 128, 256], 

„ -. * 3?min X max CFL Ax 

aw = 0, x max = 1, Ax = , CFL = 0.5, At = . , 

N C eii f 1/2 C 

IC for Gaussian Pulse: E®, F r ° = exp [ — (u(x — /i)) 2 ) , v = 20, \i = 0.3, 

1 0.2 < x < 0.4 
otherwise 

Since the Gaussian pulse results from smooth initial data, one expects R\ = 2.0. However, 
the square wave results from discontinuous initial data and one expects R\ ~ 0.67. This 
claim is true for all second order spatially accurate numerical methods when applied to an 
advection-type problem (u t + au x = 0) [I]. 



IC for Square Pulse: £° F r ° 



4.3 Weak Equilibrium Diffusion Limit 

In the weak equilibrium diffusion limit, r > 0(1) and the radiation subsystem reduces to 
Equations [HI O The optical depth suggests the range of total opacities for which diffusion is 
observed: if r = cr t £diff > 1, then one expects diffusive behavior for a t > 1/^diff- Additionally, 
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Figure 2: Gaussian pulse in free streaming 
limit, t = 4 x 10- 6 = 0.4 (x max - x min )/C. 
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Figure 3: Square pulse in free streaming 
limit, t = 4 x 10- 6 = 0.4 (x max - x min )/C. 
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1.5 
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8.6E-4 


2.1 


3.1E-2 


1.4 


8.6E-4 


2.1 


3.1E-2 
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Table 2: Errors and convergence rates for Gaussian pulse in free streaming limit. Errors 
were obtained through analytic comparison. t = 4x 10~ 6 = 0.4 (rc max — x m i n )/C. 



N cM L^Er) Rate ^lC^r) Rate 

32 6.0E-2 - 6.0E-2 

64 4.2E-2 0.5 4.2E-2 0.5 

128 2.6E-2 0.7 2.6E-2 0.7 

256 1.5E-2 0.8 1.5E-2 0.8 



Table 3: Errors and convergence rates for square pulse in free streaming limit. Errors were 
obtained through analytic comparison, t — 4 x 10 -6 = 0.4 (x max — x min )/C 



14 



M. Sekora, J. Stone 




Equations fTTl [T2l set the time scale t^is and length scale idis for diffusion, where tdm ~ £&\sl D 
and .D = fC/a t for the radiation subsystem. Given a diffusion problem for a Gaussian pulse 
defined over the entire real line (ut — Du xx = 0), the analytic solution is given by the method 
of Green's functions: 

u(x, t) = [ f(x)G(x,t]x,0)dx=—— — , exp ^ 



(ADtu 2 + 1)1/2 ^ ADtu 2 + 1 



Parameters: 



C = 10 5 , a a = 40, a t = 40, / = 1/3, T 4 = E r , 
N ceU = [320, 640, 1280, 2560], 

2-min — 5, x max — 5, /\x — , CFL = 0.5, At 



IC for Gaussian Pulse: 



N cell ' ' ' /i/ 2 C 

E r ° = exp - /i)) 2 ) , v = 20, /i = 0.3, 



One's intuition about diffusive processes is based on considering an infinite domain. So 
to minimize boundary effects in the numerical calculation, the computational domain and 
number of grid cells were expanded by a factor of 10. In Figures HI El one observes the 
diffusive behavior expected for this parameter regime. Additionally, the numerical solution 
compares well with the analytic solution for a diffusion process defined over the entire real 
line (Equation |43|) . However, diffusive behavior is only a first order approximation to more 
complicated hyperbolic-parabolic dynamics taking place in radiation hydrodynamics as well 
as the radiation subsystem. Therefore, one needs to compare the numerical solution self- 
similarly. In Table 4, one sees first order convergence when a hyperbolic time step Ath = 
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N ce ll 


Li{E r ) 


Rate 




Rate 


Li(F r ) 


Rate 


Loo{F r ) 


Rate 


320 


8.9E-3 




4.5E-2 




1.1E-3 




3.7E-3 




640 


6.6E-3 


0.4 


3.4E-2 


0.4 


8.3E-4 


0.4 


3.1E-3 


0.2 


1280 


3.4E-3 


1.0 


1.6E-2 


1.1 


4.1E-4 


1.0 


1.4E-3 


1.2 


2560 


1.6E-3 


1.1 


7.1E-3 


1.1 


1.9E-4 


1.1 


6.0E-4 


1.2 



Table 4: Errors and convergence rates for E r , F r in the weak equilibrium diffusion limit. 
Time was advanced according to a hyperbolic time step: At h = C J 1 L /2 ^ X - Errors were 
obtained through self-similar comparison, t — 4 x 10~ 6 . 



N ce ll 


L^E,.) 


Rate 


Loo(E r ) 


Rate 


Li(F r ) 


Rate 


Loo(F r ) 


Rate 


320 


1.7E-2 




8.3E-2 




2.0E-3 




7.9E-3 




640 


5.0E-3 


1.7 


2.5E-2 


1.7 


6.0E-4 


1.7 


2.0E-3 


2.0 


1280 


1.1E-3 


2.2 


5.1E-3 


2.3 


1.3E-4 


2.3 


3.6E-4 


2.4 


2560 


2.5E-4 


2.1 


1.2E-3 


2.1 


2.8E-5 


2.2 


7.4E-5 


2.3 



Table 5: Errors and convergence rates for E r , F r in the weak equilibrium diffusion limit. 
Time was advanced according to a parabolic time step: At p = CFL 2 ^ X ^ ■ Errors were 
obtained through self-similar comparison, t = 4 x 10 -6 . 
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Nceii L±(E r ) Rate L 00 (£ , r ) Rate 

320 2.2E-3 - 1.8E-2 

640 5.3E-4 2.1 5.6E-3 1.6 

1280 1.3E-4 2.0 1.5E-3 1.9 

2560 3.3E-5 2.0 3.8E-4 2.0 



Table 6: Errors and convergence rates for E r in the strong equilibrium diffusion limit. 
Errors were obtained through self-similar comparison, £ = 4 x 10~ 6 . 



C fi/2£ X is used; while in Table 5, one sees second order convergence when a parabolic time 
step At p = CFL 2 ^ X " > is used. This difference in the convergence rate results from the 
temporal accuracy in the numerical solution. In the weak equilibrium diffusion limit, the 
Godunov method reduces to a forward-time/centered- space discretization of the diffusion 
equation. Such a discretization requires a parabolic time step At ~ (Ax) 2 in order to see 
second order convergence because the truncation error of the forward-time/centered-space 
discretization of the diffusion equation is 0(At, (Ax) 2 ). 

4.4 Strong Equilibrium Diffusion Limit 

In the strong equilibrium diffusion limit, r ^> 0(1). From Equations [13], [HJ F r — > for all 
time and space while E r = E®. 

Parameters: 

C = 10 5 , <7« = 10 6 , a t = 10 6 , / = 1/3, T 4 = E r , 
Ncdi = [320, 640, 1280, 2560], 

_C C A — Xmin ~ XmSLX /inr _nr A . CFL Ax 

N ceU fy*c 

{ F,° = 

IC for Gaussian Pulse: 



E° r = exp {-{u(x - fi)) 2 ) , v = 20, p. = 0.3, 

F = _±dEl = 2fu 2 {x-y) p0 



In this test, the numerical solution is held fixed at the initial distribution because a a , at are 
so large. However, if one fixed £diff and scaled time according to idiff ~ ^d?s/ D = £j* s a t / fC, 
then one would observe behavior similar to Figures HI [5j This test illustrates the robustness 
of the Godunov method to handle very stiff source terms. 
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5 Conclusions and Future Work 

This paper presents a Godunov method for the radiation subsystem of radiation hydro- 
dynamics that is second order accurate in both time and space, unsplit, asymptotically 
preserving, and uniformly well behaved. Moreover, the method employs familiar algorithmic 
machinery without a significant increase in computational cost. This work is the starting 
point for developing a Godunov method for full radiation hydrodynamics. The ideas in 
this paper should easily extend to the full system in one and multiple dimensions using a 
MUSCL or CTU approach [24]. A modified Godunov method that is explicit on the fastest 
hyperbolic scale (radiation flow) as well as a hybrid method that incorporates a backward 
Euler upwinding scheme for the radiation components and the modified Godunov scheme for 
the material components are under construction for full radiation hydrodynamics. A goal 
of future research is to directly compare these two methods in various limits for different 
values of cja^. Nevertheless, one expects the modified Godunov method that is explicit on 
the fastest hyperbolic scale to exhibit second order accuracy for all conservative variables 
and the hybrid method to exhibit first order accuracy in the radiation variables and second 
order accuracy in the material variables. Work is also being conducted on applying short 
characteristic and Monte Carlo methods to solve the photon transport equation and obtain 
the variable tensor Eddington factors. In the present work, these factors were taken to be 
constant in their respective limits. 
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